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FOREWORD 


The  work  reported  herein  was  performed  at  the  Indian  Head  Division,  Naval  Surface  Warfare 
Center  (NSWC).  The  first  part  of  this  paper,  which  concerns  the  derivation  of  posterior  marginal 
distributions  for  a  prior  consisting  of  a  mixture  of  ordered  Dirichlet  distributions,  was  presented  on  6 
November  1982  at  the  Virginia  Polytechnic  and  State  University  symposium,  “Reflections  on  Bayesian 
Approaches  in  Operations  Research,  Probability,  and  Statistics,”  in  Blacksburg,  VA.  At  the  symposium 
the  author  learned  that  the  posterior  marginals  for  an  ordered  Dirichlet  prior  had  been  published  a  year 
earlier  by  Damon  Disch  (1981).  In  November  1984  the  author  revised  the  original  paper  and  included 
recursive  relationships  that  enable  the  posterior  marginal  distributions  derived  earlier  to  be  calculated. 
These  results  were  not  published,  but  they  have  been  applied  in  an  interactive  computer  code  MBR 
written  by  the  author  and  Mr.  Patrick  O’Neal  at  NSWC  (White  Oak  Laboratory)  circa  1988.  This  code 
was  recently  revised  and  rewritten  in  Mathcad  and  published  as  IHTR  2323. 
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INTRODUCTION 


This  paper  concerns  the  problem  commonly  associated  with  bioassay,  but  there  are  broader  areas 
of  application  such  as  accelerated  life  testing,  sensitivity  testing,  and  military  damage  analysis.  We  will 
be  concerned  with  items  subjected  to  differing  levels  of  stress  and  with  the  problem  of  estimating  the 
probabilities  of  failure  associated  with  the  stress  levels.  Emphasis  is  placed  on  the  attainment  of  interval 
estimates.  The  approach  we  take  is  Bayesian  and  the  model  we  develop  is  an  extension  of  that  published 
by  Ramsey  (1972),  whose  work  was  brought  to  the  attention  of  the  author  by  Professor  N.  D. 
Singpurwalla  of  George  Washington  University. 

The  number  of  failures  y  at  the  zth  stress  level  S  is  taken  to  be  binomially  distributed 


v  ~  b(nj ,  Pi )  (1) 

where  «,•  is  the  number  of  tests  and  pt  is  the  unknown  (random)  probability  of  failure.  We  let  the  total 
number  of  stresses  involved  be  M.  It  is  assumed  that  the  unknown  p  values  underlying  the  tests  satisfy 
the  same  (complete)  ordering  restrictions  as  the  stresses.  These  we  are  free  to  write  as 

Sl<S2<-<SM  (2) 

and 

P\  <  Pi  <  " '  <  Pm  •  (3) 

In  the  discussion  that  follows  we  develop  a  joint  prior  for  the  p  values  that  is  consistent  with  the 
ordering  (3).  This  prior  is  related  to  that  proposed  by  Ramsey  (1972),  but  is  considerably  less  restrictive. 
Specification  of  the  prior  is  achieved  by  specifying  its  marginals,  which  can  be  accomplished  in  a 
variety  of  ways.  The  method  currently  used  by  the  author  is  to  obtain  the  user’s  judgments  as  to  the 
modes  (most  likely  values)  and  5th  and  95th  percentiles  (uncertainty  limits)  at  each  of  the  stress  values. 
Usually  these  can  be  obtained  in  the  forms  of  modal  and  limits  curves  spanning  the  stress  values  of 
interest.  The  complete  marginals  can  then  be  supplied  by  the  statistician  in  a  manner  that  is  consistent 
with  the  user’s  choices.  As  will  be  shown,  the  marginals  must  satisfy  certain  conditions  that  amount  to 
very  natural  restrictions  on  the  forms  of  the  marginal  distribution  functions. 

In  the  final  sections  we  develop  expressions  for  the  posterior  marginals.  From  these  expressions, 
the  user  can  obtain  new  values  of  the  marginal  modes  and  percentiles  as  modified  by  the  data. 
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CONSTRUCTION  OF  THE  PRIOR 


Ramsey’s  prior  was  a  form  of  the  ordered  M-variate  Dirichlet  distribution,  which,  following  Wilks 
(1962,  p.  182),  we  can  write  generally  as 


gj(P\<P2<---<PM)  = 


r  (aiy  -I  +aM+lj)  T-r1  , 
r(al7.)-r(aM+w)}i  Pi 


-Pi- 1)  ]  » 


(4) 


where  0  =  pG  <  py  <  p2  <•■■  <  Pm  <  Pm+ 1  =1  and  ot,y  >0  for  all  ij.  Our  prior  will  consist  of  a 
mixture,  or  convex  combination,  of  (4),  viz., 


J 

g(p)  =  ^gj(p), 

.7=1 


(5) 


where  <j)y  >  0 ,j=\ ,  2,  ...,  J,  and  ^())y  =1.1 

j 

Now,  we  want  to  assign  values  to  the  parameters  {a;y }  in  (4)  and  (5)  by  choosing  the  shapes  of  the 
M  marginals  of  g(p)  .  These  we  require  to  be  of  the  fonn 


J 

g(Pi )  =  X't’./  SjiPi) »  z'  =  1.2, ...,M, 
.7=1 


where 


gj(Pi)  = 


r  (ciij  +  by ) 


r  (al7)r  (6./) 


(1-P/A"1 


(6) 


(7) 


Hence,  our  joint  prior  consists  of  a  mixture  of  ordered  M-variate  Dirichlet  distributions  with  marginals 
that  are  mixtures  of  betas.  So  that  (6)  and  (7)  are  bounded,  we  will  require 

djj  >  1  and  kj  >  1  •  (8) 


1  this  choice  was  motivated  by  Mazzuchi’s  success  (Mazzuchi  and  Singpnrwalla.  1981)  in  representing  the  moments  of 
the  posterior  marginals  for  (4),  i.e.,  Ramsey’s  prior,  and  the  need  for  a  richer  class  of  priors. 
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For  the  above  forms  to  hold,  Wilks  (1962,  theorem  7.7.6)  gives  the  following  conditions 


£=\ 


aU=La(j 


and 


M+ 1 

bij  =  Z a 


£=i+ 1 


tj 


(9) 


(10) 


for  i  =  1,2,...,  M  and  j  =  1,2,...,  J.  From  (9)  we  find 


a  1/  a  1  j 
a2  j  ~  a2j  ~  a\  j 


(11) 


aMj  -aMj~aM-\,j 


Hence  a1y,a2y,...,ajWy  are  fully  detennined  by  aij,a2j .  And  from  (10)  we  find 


aM+\,j  =bMj  (12) 

and  by  =bMj  +  aMj-aij  ,  i  =  (13) 

The  latter  result  is  found  by  substituting  (11)  and  (12)  into  (10).  Also,  within  the  additive  constant, 
bj^j  j ,  the  b  parameters  are  determined  by  the  a  parameters.  From  these  results  we  find  that  the  beta 

densities  of  (7)  are  necessarily  of  the  form 


gj(Pi)  «  P, 


aij~  1 


(1  _p.^bMj+aMj-aij-1 


(14) 


It  is  convenient  to  reparameterize  (14)  in  terms  of  its  mode  p*j  and  precision  index  [!  ,  which  are 


expressed  as 


*  az/  1 

PU=~V~ 


(15) 


P j  =  bMj  +aMj- -2, 


(16) 


from  (8).  Substituting  these  into  (14),  we  obtain 
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(  *  ,  *  APj 

Pv  n  _ 


gj(Pi)  «  P  11  (1  -Pi) 


J 


Expression  of  the  Dirichlet  parameters  in  terms  of  [1  .  and  p yields 


a//  fijiPij  Pi-l,j)+&i,M+l 


(17) 


(18) 


for  i  =  1,2,...,  M+ 1  and  j  =  1,2...,  J.  Here,  as  above,  we  have  defined  pgy  =  0  and  Pm+IJ  =  1-  Note  that 
the  presence  of  the  Kronecker  8  ’s  in  (18)  makes  a M+l,j  anc*  a\  j  depart  from  the  values  assigned  by 
Ramsey,  which  appear  to  be  in  error  (cf.  p.  844). 

Thus,  by  the  above  construction  we  are  able  to  specify  a  suitable  form  for  the  joint  prior  if  we  are 
able  to  represent  our  marginal  priors  by  functions  of  the  form 


where  B  1  (u,  v)  =  —  — +  V  ,  and  i  =  1,2,...,  M. 

r(w)T(v) 

We  now  show  a  method  for  choosing  the  parameters  of  (19)  that  will  permit  the  representation  of  a  wide 
class  of  marginal  distributions. 


Assignment  Of  Parameters 

The  assignment  problem  can  be  stated  as  follows.  We  wish  to  assign  values  to  the  parameters 
cj )j,  (3  j  ,j  =  1,2,...,  J  and  p*\ , /?*2 ,. . .,  p*j  ,  i  =  1,2,...,  M  in  such  a  manner  that  the  constraints  <j)y  >  0, 

(3y  >  0,  a ij  =  $j(p*j  ~ P*i-\ ,j)+&i,M+l  +8/1  >  0,j  =  1,2,...,  J,  i=  1,2,...,  M+l  and  ^cj)y  =1  are  satisfied 
and  g(Pi),  i  =  1,2  ,...,M,  are  close  representations  of  the  actual  prior  marginals.  In  the  following  let  us 

denote  the  densities  and  distribution  functions  of  the  actual  prior  marginals  by  g*(-)  and  G*(-), 

respectively.  It  is  assumed  that  such  functions  are  available  upon  consultation  with  the  user  as  described 
earlier. 
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It  is  noted  that  choices  of  the  and  [}  /  parameters  affect  all  of  the  marginals  represented  by  (19) 
simultaneously,  but  that  sets  of  p*j  parameters  may  be  independently  assigned.  This  and  the  recognition 
that  distributions  may  be  regarded  as  describing  the  concentrations  or  densities  of  units  of  probability 
suggests  the  following  assignment  plan.  For  j  =  1,2,  . . .,  /  let 

<(>.,•=  1//  (20) 

and  (3y  =  (3  (21) 

where  p  is  a  large  number  (e.g.,  500  or  1000).  Then  for  the  /'th  marginal  we  choose  the  p*j  ,j  =  1,2,  ..., 
/as  percentiles  of  g* (■)  as  follows: 


G’(p’)  =  j/(J  + 1).  (22) 

This  process  is  repeated  for  all  marginals  g* (•) ,  i  =  1,2,...,  M.  The  interval  between  adjacent  />*•  values 
corresponds  to  equal  and  constant  units  of  probability. 

This  method  of  approximating  the  desired  distribution  shapes  is  very  similar  to  that  used  in  pattern 
recognition  theory  and  found  in  the  theory  of  Parzen  estimators  (see,  e.g.,  Fukunaga,  1972,  p.  166).  By 
increasing  the  value  of  /,  the  accuracy  of  the  approximation  can  be  increased  arbitrarily.  The  size  of 
the  P  parameter,  which  controls  the  width  of  the  beta  “kernels”  (i.e.,  components  of  the  sum),  should  be 
chosen  so  that  all  kernels  overlap  to  some  extent  (see  Meisel,  1972,  p.  107).  In  the  author’s  experience 
reasonable  choices  for /may  be  less  than  50. 

A  review  of  the  constraint  requirements  stated  at  the  beginning  of  this  section  shows  that  the  only 
constraints  requiring  our  attention  under  this  plan  are  those  on  the  a;/-  ’s,  namely 

aU  =  P  j(Pij  ~  Pi-l,j)  +  ^i,M+l  +5/i  >  0  ,  i  =  1,2,  ...,M  + 1 ,  j  =  1,2, ...,/.  It  is  of  interest  to  see  how  these 
affect  the  shapes  of  the  marginals  that  can  be  represented  by  (19). 

An  obvious  set  of  sufficient  constraint  conditions  is  given  by 


P*-\j  <  p*j  ,  i  =  1,2,  ...,M  +  1,  y'  =  l,2,...,/. 


(23) 


Now,  proceeding  pairwise,  the  plan  requires  that  G*_\{p*_\  j)  =  G*  (p*  )  =  j  /(J  =  1)  ,  which  is  illustrated 
in  Figure  1.  Hence,  the  constraint  p*-\j  <  p*j  implies 
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G*-\  (p)  >  G* (p) ,  and  p  e  .  (24) 

This  may  be  equivalently  stated  as 

Pr(Pi-\  >  P)  <  MPi  ~P)  ,  P^  [Pi-l,j ,  Pij }  •  (25) 


Figure  1.  Illustration  of  Modes  Assignments 


By  considering  other  values  of  j  while  allowing/,  P  — >  oo  ,  we  can  extend  the  region  of  validity  of  (24) 
and  (25)  to  the  entire  unit  interval  (0,1).  Thus,  the  full  set  of  constraints  given  in  (23)  suggests  the 
following  restrictions  on  the  marginals  that  can  be  approximated  by  (19): 


G\(P)  >  Gl{p )  >  •  •  •  >  G*m  (p)  for  all  p  €  (0,1)  , 


(26) 


which  is  equivalent  to 


Pr( P]  >  p)  <  Pr(/»2  >  p)  <  •  •  •  <  Pr {pM  >  p),  for  all  p  e  (0,1).  (27) 

Both  relations  (26)  and  (27)  express  conditions  that  must  be  satisfied  by  the  prior  marginals.  The  latter, 
which  are  referred  to  as  the  conditions  for  stochastic  ordering,  have  an  intuitively  meaningful 
interpretation. 
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Other  Features  Of  The  Prior 

Interesting  properties  of  Ramsey’s  prior  that  were  listed  in  his  paper  (1972)  appear  to  carry  over  to 
the  mixed  prior.  The  conditional  distribution  of  pt  given  pi_\  and  pj+\  is  a  mixture  of  translated  betas  on 

the  interval  (Pi-\,Pi+\)  ■  And  in  the  limit  as  J  — »  oo  and  p  — »  qo,  g(p)  can  be  written  as 

g(p)  =  g(P\)  g(P2  \P\)g(P3\Pl)---g(PM  \Pm-\)>  (28) 

where  the  conditionals  are  sums  of  translated  betas  respectively  over  •••,  (Pm- i>1)  (also 

see  Kraft  and  Van  Eeden,  1964). 
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THE  POSTERIOR  DISTRIBUTION 


Since  the  data  are  assumed  to  be  binomially  distributed  and  acquired  independently,  the  likelihood 
function  is 


M 

L{p  I  y)  =  n 


Pp  (!-  Pi) 


7=t  w 


>y  -y, 


Hence,  the  joint  posterior  density  is  written  as 

f(p\y)ccL(p\y)  g(p)  , 


where  g(p)  is  obtained  from  (6).  Substituting  this,  we  get 


(29) 


(30) 


J 

r  m  i 

a  -pm  n  p>‘  a  -  Pi  (p;  -  n- a 

f(p  1  y) «  Yj$jKj 

,/=i 

i=i 

(31) 


where 


_  T  (a-ij  +---+aM+1j) 

J  r(a1_/)---r(a  M+ij) 


(32) 


The  marginals  of  the  posterior  distribution  are  thus  obtained  by  perfonning  the  integrations  indicated  in 


Pi  P  2  1  1 

$jKj  J  dPi- 1  dPi+ 1  ' ' '  J  dPM  {'} 


Pm -i 


(33) 


where  {•}  is  the  bracketed  term  in  (3 1)  above. 

Integration  of  (33)  can  be  achieved  by  expanding  the  various  binomial  terms  in  a  systematic 
fashion.  The  process  has  been  described  by  Disch  (1981)  in  considerable  detail  (for  J  =  1).  Using  a 
notation  similar  to  that  of  Disch,  the  result,  for  i  =  1,2,  . . .,  M,  can  be  expressed  as 
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J  "1  «i-l 


M  fn  ^ 


f^P\y}^Y^jKj  Z  "'  Z  (-D^n  0  B  l(£rj+K>ar+ \,j) 


j= 1  <i=0  <m=0 


,=1Wry 


A+l  >;M 


z  -  z  (-dAm  n 

fi+1=0  ^M=0  5=1+1  V"'‘?y 


( 


(34) 


where 


qi=l-Pi,  nk=nk- yk 


^ \rj  ~  (kA:  —  '/c  ’  ^0  —  0 

A:=l  A:=l 


M  M 

fi  rj  =  (w£  “*"aA:+l,j/ )  ’  As  =  ^  k  ■>  ^~M  +1  =  0 

k=s  k=s 


and  the  B  Unction  was  defined  earlier  in  connection  with  (19).  Clearly,  from  (34)  one  can  easily 
obtain  forms  for  the  M  posterior  marginal  distribution  functions  in  terms  of  incomplete  beta  functions. 

Interpolation  and  Extrapolation 

Expression  (34)  is  readily  extended  to  non-observational  stresses  (i.e.,  stresses  for  which  no  data 
exist)  by  setting  the  associated  values  of  h,  and  v,  to  zero.  In  this  way  we  can  interpolate  or  extrapolate 

the  forms  of  the  posterior  marginals  at  stress  values  other  than  those  at  which  data  were  collected 
provided  the  prior  marginals  at  these  points  have  been  specified. 

Computation  of  the  Posterior  Density 

In  this  section  we  develop  a  recursive  procedure  by  which  (34)  can  be  rapidly  calculated. 
Disch  (1981)  explored  approximate  methods  of  computation  after  pointing  out  the  profound  inefficiency 
of  straightforward,  brute-force  approaches.  Antoniak  (1974)  reported  similar  difficulties. 

Equation  (34)  may  be  viewed  as  the  sum,  over  j,  of  the  product  of  two  sums,  the  first  being  a  sum 
over  all  (z-l)-tuples  (t\,£ 2,—,f/-i)and  the  second  a  sum  over  all (M - z)-tuples  (f  j+i,f/+2,  —  ,^M)  •  A 
useful  alternative  representation  involves  the  summations  in  which  the  indices  l\,i and 
are  constrained.  We  then  write  (34)  as 
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J 


f(Pi\y)azYu$jK'j 

7=1 


— < 

n 

i- 1 


Z  (-1)  C p(i,  j , k) pj 
£=0 


£//  +^-l 


z  (-1 


k'= 0 


(35) 


where 


Cp(i,j,k)  =  Z  -  2  f] 


A 


£1=0  7(_,  =0  r=l  V^rJ(Cry) 


A 


subject  to  the  constraint  =  k  ,  and 


.>7+1  >"m  M 


Cq  (A  j  ■,  k’)  =  Z  Z  II 


ys 

L 


'a. 

A 

V 

'a. 

Vl=0  V=0  ,S'=I  V  W 

J 

subject  to  the  constraint  A;+i  =  k' .  Here,  we  have  made  use  of  Pochhammer’s  symbol 

{z)k  =T(z  +  A:)/r(z)  =  z(z  +  l)---(z  +  A:-l),  (z)Q=l  , 
where  k  is  an  integer,  and  we  have  defined 


(36) 


(37) 


C/ j  —  Aj  +ar+lJ  5  ^77  —  fis/  ~1~ ^ ,sy  ■ 
7-1  M 

V,  =  IX  ,  V,  =  E  A, 


?'=i 


and 


K  ■  =  K  ■  TT  r  ^/y  ^r  (-CXr+17^  TT  r  (l1.v/)r  (asj) 

J  r= 1  r  (^77  +  ar+ly)  s,=f+|  f  Olsy  +  a,s/) 


As  a  consequence  of  (36)  and  (37),  we  find  Cp  (1,  j,  0)  =  Cq  (M,  j,  0)  =  1 . 

The  simplification  of  equations  (36)  and  (37)  can  be  effected  by  making  index  transformations 
which  pennit  the  factoring  of  the  summands  through  the  summations.  We  first  examine  equation  (36). 
Consider  the  1  -to- 1  and  onto  transfonnation  of  summation  indices  from  the  set  to  the  set 
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A,i,A,2,...,A,/— i  with  transformation  relations  ir  =Xr- Zr_]  ,  r  =  1,2,...,/  —  1,  =  0 .  Under  these 

transformations  (36)  can  be  expressed  as 

min (k,nr_2  )  min ( V2  >»r_3  ) 

Cp(i,j,k)=  Yj  Y 

hj_2  =max(£-7?!_1 ,0)  Xj_  3  =max(A,_2  -n,_2  >0) 

min(),2,«f )  /-I  / 

2  n  I, 

^1=max(^2“«2’0)  r=l  x  7  7  ^ 


where  A,;_i  =  k  and  the  summation  limits  arise  from  the  constraints 


~)-r  —  nr  <  Xr_i  <  ~>-r 


\  >0 


kr  <  tiy  . 


By  factoring  out  the  product  terms,  we  obtain 

fe-h/)/ 


min (k,n~  ) 

i-2 


Cp(i,j,k)  = 


Z 


k  \-2  =max(&-ni_1 ,0) 


^  »/-l  A 

\k-\_2J 


min()w_2,«^3 ) 


z 


'  »/-2  A 


(^'-2,y)x._2  U-3  =rnaxiU-2  ~ni-2  >o)  2 


X  •••  X 


(^2’A2 

min(Z2,7?~) 

V 

'  «2  " 

/«i  N 

2—i 

>4=max(Z2-/?2,0) 

vz2  -Zly 

(<*4 

v^U 

(39) 


By  inspection  of  (39),  we  find  the  following  recursive  relationship  between  the  Cp  coefficients: 
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c p  0?  j ? 


min(A:,/7<  ) 

_ i-2 

i 

r=max(k-ni_]  ,0) 


A »/-l A 
\k~r  J 


Cp(i-\,j,r) 


A:  =  0,1, 


h<  ,  z  =  3,  ...,M 

i- 1 


Cp(2J,k) 


M* 

A:  =  0,1, 


(40) 


Cp(hj,k)  =  1  ,  *  =  0 


Recursive  relationships  for  the  Cq  coefficients  follow  in  a  similar  manner  upon  transformation  from  the 
indices  to  the  set  A;+1,  A;+2,...,  .  We  obtain 


Cg(i,j,k)  = 


mi”^)  r*+o 


k  s=max(k-yi+\  ,0) 


Cq(i  +  \,j,s )  ,k  =  0,1,...,  v  , /=1,2,..., M-2 

1  i+ 1 


Cq(M  -  \,j,k)  = 


r  yM^ 

(v«.4i 

l  k  J 

,  k  =  0,1,...,  v 


M 


(41) 


CJM,j,k )  =  1  ,  £  =  0 


We  note  from  (9)  that  c,y  can  be  expressed  as 

i 

^>ij  =  (A/c  =  aij  +  V/ 

it=l 

where  V/“  has  been  defined  (in  a  manner  analogous  to  that  of  nf)  as  the  sum  of  all  upvalues  having  k 
less  than  or  equal  to  i.  Similarly,  from  (10)  we  note  that 

M 

fi ij  =  Yj  ("k  +^k+l,j)  =  bij  +  » 

k=i 
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where  is  the  sum  of  all  ri/(  having  k  greater  than  or  equal  to  i,  analogously  to  y,-  .  Substituting  these 
expressions  into  (35)  and  rearranging  yields 


_< 

n 

i — 1 


i+l 


f(Pi  I  y)  «  X  (—\yC Cp(i,  j\k)  2  (-If  Cq(ijX)Pi 

7=1  A:=0  A:'=0 


ciy+yf  +k— 1  bjj+rtj  +k'~  1 


(42) 


The  constant  of  proportionality,  if  desired,  can  be  obtained  from  the  condition 


f(Pi  I  y¥Pi  =  1  , 

■00  — 


which  holds  for  any  value  of  i. 

Equation  (42)  is  a  particularly  useful  computational  fonn  for  the  posterior  marginals.  It  is  cast  in 
tenns  of  the  marginal  beta  coefficients  { <7;/- ,  /);/ }  rather  than  the  joint  distribution  coefficients  {a;/-  j  . 

This  avoids  possible  round  off  problems  associated  with  the  calculation  of  in  equation  (18) 
involving  p*j  differences. 

A  computer  program  that  uses  a  closely  related  formulation  and  contains  an  example  was 
recently  published  as  McDonald  (2003). 
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